{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mixing system\n",
    "\n",
    "Problem statement: The figure below shows a set of well-mixed mixing tanks. All the streams contain a binary mixture of substance X and substance Y. Steams A, B and C are fed into the system from an upstream process.\n",
    "\n",
    "Tanks 1 and 2 are drained by the force of gravity (assume flow is proportional to level), while the pump attached to the tank 3 output is sized such that the level in tank 3 does not affect the flowrate through the pump.\n",
    "\n",
    "You may assume that the valves in lines G, H, J and L can manipulate those flows directly.\n",
    "\n",
    "The density of substance X is ρX = 1000 kg/m3 and the density of substance Y is ρY = 800 kg/m3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img width=\"50%\" src=\"../../assets/mixing_tanks.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Steady state calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the steady state flow rates and compositions of all the streams given that\n",
    "3\n",
    "* Stream A is 1m3/h of substance X\n",
    "* Streams B and C are both 1m3/h of substance Y.\n",
    "* H=G,H=2J,J=2L."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "ρx = 1000  # kg/m3\n",
    "ρy = 800  # kg/m3\n",
    "\n",
    "A = 1*ρx\n",
    "B = 1*ρy\n",
    "C = 1*ρy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = A + B + C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = G\n",
    "J = H/2\n",
    "L = J/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = G + H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = A + L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = J - L\n",
    "I = H - J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = B + D + K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000,\n",
       " 800,\n",
       " 800,\n",
       " 1650.0,\n",
       " 3100.0,\n",
       " 5200,\n",
       " 2600,\n",
       " 2600,\n",
       " 1300.0,\n",
       " 1300.0,\n",
       " 650.0,\n",
       " 650.0)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A, B, C, D, E, F, G, H, I, J, K, L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compositions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "xA = 1\n",
    "xB = 0\n",
    "xC = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "xG = (xA*A + xB*B + xC*C)/G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "x3 = xF = xH = xI = xJ = xK = xL = xG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = xD = (xA*A + xL*L)/D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "x2 = xE = (xB*B + xD*D + xK*K)/E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Design\n",
    "Assuming all three tanks are of constant cross-sectional area of 3m2, find out what the proportionality constants should be for tank 1 and 2 so that the steady state levels will be 1 m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "A1 = A2 = A3 = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "h1 = h2 = h3 = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "k1 = D/h1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "k2 = E/h2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1650.0, 3100.0)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k1, k2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dynamic simulation\n",
    "\n",
    "Now that you have all the parameters in your system, simulate the response of the system to a sudden increase in flow rate of A from 1 m3/h to 1.5 m3/h at time 0. You should start your simulation at steady state.\n",
    "\n",
    "Assume that the level in tank 3 is also 1 m at the initial conditions. Note that the steady state relationships between H, G, J and L will not hold over the whole simulation. Simply set them to their steady state values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our states will be the total masses and mass of X in each tank. Let's find the initial values at steady state first:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find volumes via tank geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "V1 = A1*h1\n",
    "V2 = A2*h2\n",
    "V3 = A3*h3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Masses from volumes - assume ideal mixing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "M1 = V1/(x1/ρx + (1 - x1)/ρy)\n",
    "M2 = V2/(x2/ρx + (1 - x2)/ρy)\n",
    "M3 = V3/(x3/ρx + (1 - x3)/ρy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = [M1, M2, M3, M1*x1, M2*x2, M3*x3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dMdt(t, y):\n",
    "    M1, M2, M3, M1x1, M2x2, M3x3 = y\n",
    "\n",
    "    if t <= 0:\n",
    "        A = 1*ρx\n",
    "    else:\n",
    "        A = 1.5*ρx\n",
    "    \n",
    "    xD = x1 = M1x1/M1\n",
    "    xE = x2 = M2x2/M2\n",
    "    xF = x3 = M3x3/M3\n",
    "    \n",
    "    V1 = M1*(x1/ρx + (1 - x1)/ρy)\n",
    "    V2 = M2*(x2/ρx + (1 - x2)/ρy)\n",
    "    V3 = M3*(x3/ρx + (1 - x3)/ρy)\n",
    "\n",
    "    h1 = V1/A1\n",
    "    h2 = V2/A2\n",
    "    h3 = V3/A3\n",
    "\n",
    "    xH = xI = xJ = xK = xL = xG = x3\n",
    "    \n",
    "    D = k1*h1\n",
    "    E = k2*h2\n",
    "\n",
    "    dM1dt = A + L - D\n",
    "    dM2dt = B + D + K - E\n",
    "    dM3dt = C + E + I - F\n",
    "\n",
    "    dM1x1dt = xA*A + xL*L - xD*D\n",
    "    dM2x2dt = xB*B + xD*D + xK*K - xE*E\n",
    "    dM3x3dt = xC*C + xE*E + xI*I - xF*F\n",
    "\n",
    "    return dM1dt, dM2dt, dM3dt, dM1x1dt, dM2x2dt, dM3x3dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We expect t=0 to give zero derivatives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0,\n",
       " 4.547473508864641e-13,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 4.547473508864641e-13,\n",
       " -4.547473508864641e-13)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dMdt(0, y0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And for other times to give non-zero derivatives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500.0,\n",
       " 4.547473508864641e-13,\n",
       " 0.0,\n",
       " 500.0,\n",
       " 4.547473508864641e-13,\n",
       " -4.547473508864641e-13)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dMdt(1, y0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = scipy.integrate.solve_ivp(dMdt, (0, 10), y0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "  message: 'The solver successfully reached the end of the integration interval.'\n",
       "     nfev: 68\n",
       "     njev: 0\n",
       "      nlu: 0\n",
       "      sol: None\n",
       "   status: 0\n",
       "  success: True\n",
       "        t: array([0.00000000e+00, 1.00000000e-04, 9.32930762e-04, 9.26223838e-03,\n",
       "       9.25553145e-02, 9.25486076e-01, 1.86268269e+00, 3.13373601e+00,\n",
       "       4.83441945e+00, 7.23515873e+00, 9.99326358e+00, 1.00000000e+01])\n",
       " t_events: None\n",
       "        y: array([[2828.57142857, 2828.61687016, 2829.03321934, 2833.18623756,\n",
       "        2873.68652849, 3191.18414035, 3411.4077485 , 3576.81492267,\n",
       "        3679.20042145, 3732.50482985, 3751.73742159, 3751.76206113],\n",
       "       [2657.14285714, 2657.14285827, 2657.14297442, 2657.1545728 ,\n",
       "        2658.2647728 , 2730.90792703, 2850.15605228, 2978.77726282,\n",
       "        3075.5278406 , 3130.61308859, 3150.4898384 , 3150.51308613],\n",
       "       [2600.        , 2600.        , 2600.00000004, 2600.00003725,\n",
       "        2600.0360844 , 2626.36069917, 2755.48727041, 3096.98554741,\n",
       "        3748.1911919 , 4840.17117488, 6180.11425868, 6183.43458116],\n",
       "       [2142.85714286, 2142.90258442, 2143.31893191, 2147.4717797 ,\n",
       "        2187.95586832, 2505.19417937, 2729.14975288, 2908.4244449 ,\n",
       "        3034.24936155, 3113.23792305, 3148.959073  , 3149.01344238],\n",
       "       [1285.71428571, 1285.71428686, 1285.71440471, 1285.72617038,\n",
       "        1286.84986914, 1359.78761856, 1480.91665298, 1618.72961343,\n",
       "        1733.52304584, 1810.12997441, 1844.51567419, 1844.56643364],\n",
       "       [1000.        , 1000.        , 1000.00000004, 1000.00004033,\n",
       "        1000.03821936, 1022.64243623, 1116.4204006 , 1327.1016063 ,\n",
       "        1689.48149947, 2263.4913374 , 2943.38602418, 2945.04374286]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the composition of stream G as well as the compositions and levels in all three tanks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "M1, M2, M3, M1x1, M2x2, M3x3 = sol.y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = M1x1/M1\n",
    "x2 = M2x2/M2\n",
    "x3 = M3x3/M3\n",
    "\n",
    "V1 = M1*(x1/ρx + (1 - x1)/ρy)\n",
    "V2 = M2*(x2/ρx + (1 - x2)/ρy)\n",
    "V3 = M3*(x3/ρx + (1 - x3)/ρy)\n",
    "\n",
    "h1 = V1/A1\n",
    "h2 = V2/A2\n",
    "h3 = V3/A3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1519b93c50>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxV5bn//c+VeYIkkCBIGMKkKBXFiFgnrBNoBfrTttqjx05y9Ger7Xk8rcM59RyLtR5bW/tURKoU6tjzWLV4xAnbiiMySFBEDQSQMJiQQCCQeV/PH3sHNiEhG9jJTjbf9+u1X3uve91rrWuZ+M1irbXXbe6OiIjEr4RYFyAiIp1LQS8iEucU9CIicU5BLyIS5xT0IiJxTkEvIhLnOgx6MxtkZn83s9VmtsrMbm6jj5nZ78xsjZmtNLNxYfOuNbOS0OvaaO+AiIgcnHV0H72ZDQAGuPtyM+sFLAOmufvHYX0uAX4IXAKcDjzg7qebWR9gKVAEeGjZU919e6fsjYiIHKDDI3p33+Luy0OfdwGrgYGtuk0F/uRB7wE5oT8QFwOvuXtVKNxfAyZFdQ9EROSgkg6ls5kNBU4BFreaNRDYGDZdFmprr72tdU8HpgNkZmaeevzxxx9KaSISxxoDjWyp2cKuxl1kJGVwbNaxpCamxrqsbmXZsmXb3D2/rXkRB72ZZQF/AX7k7jtbz25jET9I+4GN7rOB2QBFRUW+dOnSSEsTkTgV8AD/8+n/8Jtlv2EoQ7l53M1cdfxVJJjuI2nNzDa0Ny+ioDezZIIh/4S7P9tGlzJgUNh0AbA51D6xVfs/ItmmiBzd1lev58537mR5+XLOGHAGd375TgZmtXlCQDrQYdCbmQGPAqvd/f52us0HfmBmTxO8GFvt7lvM7BXgF2aWG+p3EXBbFOoWkTjVFGhi3qp5zFwxk9SkVH5+5s+ZOnwqwSiSwxHJEf2ZwDXAh2a2ItR2OzAYwN1nAQsI3nGzBtgDfCc0r8rMfg4sCS13l7tXRa98EYknn1Z9yn+8/R+srlrNBYMv4I4Jd5CXnhfrsnq8DoPe3d+i7XPt4X0cuLGdeXOAOYdVnYgcFeqb63m4+GH++NEfyU7N5v6J93PhkAs7bXuNjY2UlZVRV1fXadvoLGlpaRQUFJCcnBzxMod0142ISLStKF/Bz975Geuq1zFl+BR+ctpPyE7N7tRtlpWV0atXL4YOHdqjTgm5O5WVlZSVlVFYWBjxcgp6EYmJPY17+N0Hv+PJ1U/SP7M/sy6YxZkDz+ySbdfV1fW4kAcwM/r27UtFRcUhLaegF5Eu987md7jr3bvYXLOZK4+/kpvH3UxmcmaX1tDTQr7F4dStoBeRLlNdX82vlv6K59c8z9DeQ5k7aS7jjhnX8YJyRBT0ItIlXt/wOjMWz2B73Xa+/6Xvc/3Y6/Xt1i6ioBeRTrWtdhu/WPwLXtvwGsf3OZ6Z589kdN/RsS7rqKLvEYtIp3B35q+dz9Tnp/LGxje4edzNPHnpkwr5MAsXLuSaa67p9O3oiF5Eom5zzWbuevcu3t78Nqf0O4X//PJ/Mix7WKzL6naKi4s55ZRTOn07OqIXkagJeICnPnmKr/31aywvX85t429j7qS5Cvl2FBcXs3XrVs4++2z69+/PwoULO2U7OqIXkahYV72O/3znP1levpwzjz2Tn53xM47NOjbWZXXov15YxcebWz+Q98iccGxv7rzsxA77FRcXM2bMGN58802effZZnnjiCS644IKo1gIKehE5Qo2BRuatmsdDKx4iLSmNGWfOYMrwKT32PvWu0tjYSFVVFbfccgsATU1N5OTkUFpayt133011dTXPPPNMVLaloBeRw7a6cjV3vnMnq6tWc+GQC7n99Nt73EPIIjny7gwff/wxY8eOJSEheAZ95cqVjBkzhmHDhvHoo49yxRVXRG1bCnoROWR7Gvcwe+Vs5q6aS05qDr+Z+BsuGBL9Uw7xrLi4mLFjx+6dXrlyJVOnTu2UbSnoRSRi7s7fPv8b9y65ly27tzBtxDRuKbql0x9CFo+Ki4sZP3783umPPvqIMWPGdMq2FPQiEpGNOzfyi/d/wVub3mJk7kjmnT1Pjy84Ar/+9a/3my4tLQWgsrKSO+64gw8++IB77rmH22478rGaFPQiclD1zfXM+XAOj3z4CEkJSfxb0b/xrdHfIilB8dEZ+vbty6xZs6K6Tv2kRKRdi8oWcc/ieyirKWPy0Mncctot9MvoF+uy5BAp6EXkAFtqtnDvknt5/fPXKcwu5A8X/YEJAybEuiw5TJEMDj4H+CpQ7u4HXCkws38D/ilsfaOB/NB4seuBXUAz0OTuRdEqXESir7G5kXkfz2P2ytkA3DzuZq494VqSEyMftk66n0iO6OcCvwf+1NZMd78PuA/AzC4DftxqAPDz3H3bEdYpIp1s8ZbF3L34btZVr+P8wefz09N+yoCsAbEuS6IgksHBF5nZ0AjXdxXw1JEUJCJdq3xPOb9a8iteWv8SBVkFPHj+g5xTcE6sy5Ioito5ejPLACYBPwhrduBVM3PgYXefHa3ticiRaQo08eTqJ5lZPJPG5kZuGHsD3x3zXdKS0mJdmkRZNC/GXga83eq0zZnuvtnM+gGvmdkn7r6orYXNbDowHWDw4MFRLEtEWlv+xXJmLJ5ByfYSzhp4FrePv51BvQfFuizpJNEM+itpddrG3TeH3svN7DlgPNBm0IeO9mcDFBUVeRTrEpGQytpK7l92P/PXzqd/Zn9+O/G3fGXwV/QAsjgXlaA3s2zgXODqsLZMIMHdd4U+XwTcFY3ticihaQ4088xnz/DABw9Q21TL98Z8j+knTScjOSPWpUkXiOT2yqeAiUCemZUBdwLJAO7e8vWtrwGvuvvusEWPAZ4LHSkkAU+6+8vRK11EIvHRto+Y8d4MVlWu4vT+p3P7hNs1EEg3sXDhQubNm8djjz3WqduJ5K6bqyLoM5fgbZjhbaXA2Lb6i0jnq66v5oHlD/DMZ8+Ql57Hf5/z30waOkmnaboRDSUoIocl4AGeK3mOy567jGdLnuXqE65m/rT5TC6crJDvZjSUoIgcsk+qPuHu9+5mRcUKTul3CnecfgfH9Tku1mV1by/dCls/jO46+38JJv+yw24aSlBEIrarYRcPrniQpz55ipzUHH5+5s+ZMnwKCaZ/tHdX7Q0l+Pzzz/Piiy9SXl7OjTfeyEUXXXTE21LQi/Rg7s6L617k10t/TWVtJd847hv88JQfaiCQQxHBkXdnaG8owWnTpjFt2jS2b9/OLbfcoqAXOZqt3bGWuxffzZKtSxjTdwy//8rvOTEvNuOfyqHraCjBGTNmcOONN0ZlWwp6kR5mT+MeZq2cxWOrHiMjOYP/mPAfXD7ychITEmNdmhyC9oYSdHduvfVWJk+ezLhx0RnBS0Ev0kO4Ows/X8i979/LF3u+4GsjvsaPTv0RfdL6xLo0OQztDSX4u9/9joULF1JdXc2aNWu4/vrrj3hbCnqRHmDDzg3cs/ge3t78NqNyR/Grc3/Fyf1OjnVZ0gluuukmbrrppqiuU0Ev0o3VNdXxyIePMOejOaQkpvDT037KlcdfqfFa5ZDot0Wkm3pj4xvc8/49bKrZxCWFl3BL0S3kZ+THuizpgRT0It3MpppN/PL9X/KPjf9gWPYwHr3oUcYPGN/xgiLtUNCLdBMNzQ3MXTWXP6z8A2bGj0/9MdeMvkbjtcoRU9CLdAPvbn6XXyz+Bet3rufCIRfyk9N+Qv/M/rEuS+KEgl4khr7Y/QX3Lb2PV9a/wuBeg3nogoc4a+BZsS5L4oyCXiQGGgONwfFaV8yk2Zu58eQb+c6Y75CamBrr0iQOKehFutjSrUu5e/HdrNmxhnMKzuHW8bcyqJfGa5XOo6AX6SLbardx/9L7eaH0BY7NPJYHznuA8wadp2fES6fTM0xFOllzoJknVz/JlOem8NL6l7juS9fx/LTnNSi3sHDhQq655ppO346O6EU6UXFFMXe/dzerq1YzYcAEbj/9dgqzC2NdlnQTXTWUYCSDg88BvgqUu/uYNuZPBP4KrAs1Pevud4XmTQIeABKBR9w9Ng9+FuliO+p28Nvlv+UvJX+hX3o/7jv3Pi4ecrGO4Luhe9+/l0+qPonqOo/vczw/Hf/TDvsVFxfTv39/zj77bEpKSnj88cdjNsLUXOD3wJ8O0udNd/9qeIOZJQIPAhcCZcASM5vv7h8fZq0i3V7LeK2/Wf4bahpquPaEa7nh5BvITM6MdWnSDXWboQTdfZGZDT2MdY8H1rh7KYCZPQ1MBRT0EpdWV65mxuIZrKxYybh+47hjwh2Myh0V67KkA5EceXeG9oYSXL16NQ888ADbtm3j/PPP54YbbjjibUXrHP0ZZlYMbAZucfdVwEBgY1ifMuD09lZgZtOB6QCDBw+OUlkinW9nw05+/8Hv+fOnfyYnNYe7z7qby4ZdptM0clDtDSU4evRoZs2aRSAQ4LrrrovKtqJx181yYIi7jwX+X+D5UHtbv+Xe3krcfba7F7l7UX6+ntAn3Z+788LaF7jsucv486d/5hujvsELX3uBKcOnKOSlQ20NJXjSSScBMH/+fM466yzOP//8qGzriIPe3Xe6e03o8wIg2czyCB7Bh38LpIDgEb9Ij1eyvYTvvPIdbn/rdgZmDeSpS5/ijgl30Duld6xLkx6iuLh4b7DDvqEEAaZMmcI777zDE088EZVtHfGpGzPrD3zh7m5m4wn+8agEdgAjzawQ2ARcCXzrSLcnEku7G3fz0IqHeHz142SlZHHnGXfyf0b+HxJMX0mRQ9PeUIL/+Mc/ePbZZ6mvr+eSSy6JyrYiub3yKWAikGdmZcCdQDKAu88CrgBuMLMmoBa40t0daDKzHwCvELy9ck7o3L1Ij+PuvLLhFe5bch/le8q5fOTl3DzuZnLTcmNdmsSZiRMnMnHixKiuM5K7bq7qYP7vCd5+2da8BcCCwytNpHtYV72Oexbfw7tb3mV0n9HcP/F+xuaP7XhBkW5C34wVaUdtUy1/WPkH/rjqj6QnpnPb+Nv45nHfJDEhMdaliRwSBb1IK+7O65+/zn1L7mPz7s1cNuwy/rXoX8lLz4t1aSKHRUEvEuLuvFH2BjNXzGR11WqGZw9nzsVzOK3/abEuTTqBu/fI22CDl0APjYJejnruzlub3uLBFQ+yqnIVBVkFzDhzBpcOu5SkBP0vEo/S0tKorKykb9++PSrs3Z3KykrS0tIOaTn9FstRy915d/O7PFj8ICsrVjIwayB3ffkuvjr8qyQnaEDueFZQUEBZWRkVFRWxLuWQpaWlUVBQcEjLKOjlqOPuLN66mJkrZvJB+Qf0z+zPz874GdOGTyM5UQF/NEhOTqaw8Oh5XLSCXo4qS7Yu4cEVD7Lsi2X0y+jHv5/+73xt5NdISUyJdWkinUZBL0eF5V8sZ+aKmSzeupj89HxuHX8rV4y6QoNxy1FBQS9xbUX5CmaumMm7W96lb1pffnLaT/j6qK+TlnRoF7NEejIFvcSlDys+5MHiB3l709v0SevDLUW38I3jvkF6UnqsSxPpcgp6iSurKlfx0IqHeKPsDXJSc/jRuB9x1fFXkZGcEevSRGJGQS9x4ZOqT5i5YiZ/3/h3eqf05qZTbuJbo7+lIfxEUNBLD/fZ9s94aMVDLPx8Ib2Se3HjyTfyT6P/iV4pvWJdmki3oaCXHmntjrU8VPwQr6x/hazkLG4YewNXn3C1Bv4QaYOCXnqU0upSZhXP4uV1L5OelM51X7qOa0+8luzU7FiXJp3I3Qn4vvdA6HkvgVbt7o57WDv7psPf9/VxfL/1h7WHTTsty7asJ3xd+7a1tzZvp7ZWNbeuLTU5kSljj436fz8FvfQIG3ZuYFbxLBasW0BqYirfHfNdrj3x2qNu4I9AwGkMBGhsdhqbAjQ2B2hoDk03B2gItTUFgvPD5+2bHzbdHKCxqdV0WNve6Vbrb5luCvj+gUUoyAKhesNC84CwC7SEbBtBHBaah/EMrx4rLytVQS9Hn427NvJw8cP8b+n/kpyQzD+f8M98+8Rv0ze9b6dts66xmeraRmobmg8I0sMJzwPmHRDS7Qdp6/U3BTon9RITjOREIzkxgZTEBJITE0hOajUdmp+RkrT3c2KCkZBgGJBgRoKBmWEWnN7bngAQnJ8QPt/AWtpD6zHbv1/LdHg/9tteS9+W7bVeV8v2COsTrMnYv5aWdbZba3htYbWGv4evp73171vXvnrMgj+HzqCgl25pU80mZq+czV/X/JWkhCS+NfpbfHfMdyN+JnxDU4Dq2sa9r521jeyobaB6TyPVtU3Bz6H26tpGduzZ17e+KRCVfYg0PJMSE0hPTqR3WlKoT8t8C/VJICWp1XTL/KSEA8M4ad90UqIdENQpSftPtwS2xK9IxoydA3wVKHf3MW3M/yfgp6HJGuAGdy8OzVsP7AKagSZ3L4pS3RKnttRsYfaHs3m+5HkSLIFpw7/OlCFXk+DZrNrYSHXtZqr3NOwX4uEh3fLa09B80O30Sk2id3oy2enJ5GQkM6JfFtmh6eyM4HtmStJ+4ZmccGCQKjylJ4jkiH4uwTFh/9TO/HXAue6+3cwmA7OB08Pmn+fu246oSumxAgFnV31Tm0fO1aGj7J21jZTv+YI1DfOpSngzuODO06kpP5e5q7KZy0dtrjs9OXFvUPdOT2ZQnwy+1BLWYYGdvTfQU8hOT6Z3WhJJiQld+F9BJLYiGRx8kZkNPcj8d8Im3wMO7UHJ0iO5O5W7G1i3bTelFTWUVuxm047aA46sd9Y2crDTyikpu8jot4jmrHchwcnnbI5Lm8oxfQaQc0J4SO/73HIknpqksVtFIhHtc/TfA14Km3bgVTNz4GF3n93egmY2HZgOMHjw4CiXJYerrrGZ9ZW7Ka3Yzbptu1kbCvXSihp21jXt7ZeSmMCxOWnkZKSQm5HC0L6Z+4Vz6yPrZqvmuXWP89yaZ2gKNHH5iGlcd9J1DMwaGMO9FYlPUQt6MzuPYNCfFdZ8prtvNrN+wGtm9om7L2pr+dAfgdkARUVFR9ENVbEXCDhbd9YFA3xbKMhDR+qbdtTud3tb/95pDMvPZMrJx1KYl8Ww/EyG52UxMDc9onPSVXVV/PGj2Tz9ydM0Bhq5bPhlTD9pOoN6DerEPRQ5ukUl6M3sJOARYLK7V7a0u/vm0Hu5mT0HjAfaDHrpfLvqGkOnWkKnW7btO1Kvbdx38TIzJZHC/EzGDc7lilMLGJafxbC8TArzMslMPbxfme1125m7ai5PffIU9c31XFp4Kf8y9l8Y0ntItHZPRNpxxEFvZoOBZ4Fr3P2zsPZMIMHdd4U+XwTcdaTbk4Nrag5Qtr1275H52lCor9u2m/Jd9Xv7JRgU5GYwLD+TCcP6Miw/M/jKy+KY3qlRGzC5ur6aeavm8cTqJ6htqmVy4WSuH3s9hdlHzzBuIrEWye2VTwETgTwzKwPuBJIB3H0W8DOgLzAzFA4tt1EeAzwXaksCnnT3lzthH4467s72PY17L4Ku3bbvvPnnVXtobN53riUnI5lheZmcMyp/b5APz89kcN+MTr2YubNhJ499/BiPf/w4NY01TBo6ievHXs/wnOGdtk0RaZt5N/x+cVFRkS9dujTWZXQbjc0Blm/YzqKSCt4rrWJNeQ3VtY1756ckJjCkbwaFeZnB0yz5mQwPhXpuZteOhbqrYRePr36cx1Y9xq7GXVw45EKuH3s9o3JHdWkdIkcbM1vW3neV9M3Ybmpj1R7e+KyCRZ9V8M7aSmrqm0hMMMYWZHPpSQMYlpfJ8FCoD8xJj/l94TsbdvLnT/7M3FVz2dmwk68M+gr/9+T/y3F9jotpXSKioO829jQ0sbi0am+4l27bDcDAnHQuG3ss547K48sj8uidlhzjSvepa6pjUdkiFqxbwKKyRTQGGplYMJEbTr6BE/qeEOvyRCREQR8j7s6nX+xi0WcVLPpsG++vq6KhOUBacgIThvXl6glDOGdUPsPzM6N2YTQamgJNLN6ymAXrFvD656+zu3E3eel5fPO4bzJl+BRG9x0d6xJFpBUFfRfasaeBN0u2BcO9pIIvdgbvghl1TBb/fMYQzj0un9OG9iEtuXt949PdKa4o5sXSF3l1w6tU1VXRK7kXFw+9mMmFkzntmNNITOheNcshCDRDUz00N0BzY+i9ft/npoZQW6tXUwMEmgAHD7R6hR7m3jIdSZ8j7seBbQf0a2NdB/TzVu9H0q+t/q37hbVn5sGPPoz6j1hB34mamgMUl1XvPR2zsmwHAYfs9GTOGpHHuaPyOXtUHgOy02Ndaps+2/4ZC0oX8PL6l9lUs4nUxFTOLTiXS4ZdwtkDzyYlsWsv9PZI7vsH496gbAyFaejzfkEbHrJhn9sK2g7XFcEyHp2ndUaFJez/wlq1tTzj92D9rNV7W/3a6t+qX0ICWHIH/dqo+YB+1sZ+tNEPILVzhsBU0EfZlupaFn1WwRufVfBWyTZ21jWRYDB2UA4//MpIzj0un7EFOd32yYZlu8p4ef3LvFj6Imt2rCHREplw7ARuPPlGzht0HlkpWbEuMTqa6qGuGmp3BN/rQu+12/d9bqxtP2jbPMptI5wDjR3XcqgSkiAxBRKTQ++p+z4npYTaUiApNRgc+7Unh/qHLZ8Utnz4KynlwLaW5ZJC67AESEg8wmAO6yedQkF/hOoam3l/XdXecC8prwHgmN6pTBrTn3NG5XPWiDxyMrrv0e+22m28uv5VFqxbQHFFMQCn9DuFO06/gwuHXNipg3wctkAAGna1E9Zh0+3Na6o9+PqT0iAls+1wbAm55N77h2y74Ri2TJtB206gtrcunSaTQ6SgP0zLNlTxyJvr+Nsn5dQ3BUhJTGB8YR++UTSIc0blM+qYrG51EbW1moYaXv/8dRasW8DiLYtp9mZG5Y7i5nE3M7lwctc8XKypoZ1APkhAt8yr39nBKQeDtOzgKz0n+J43av/ptJzga7/p0DLJaZ2//yJdREF/CJoDzmsff8HsRWtZ/vkOstOTuWr8YM49Lp8JhX1JT+neR1r1zfW8WfYmC9Yt4I2Nb9AQaGBg1kC+O+a7TC6czMjckYe34sZaqN4ENVvbCeg2wrquGhr3HHy9SWn7wjc9B7L6tRPWrabTcyClF6Hx60SOegr6CNQ2NPPM8jIefbOU9ZV7GNQnnf+aciJfLyogI6V7/ydsCjTx/tb3eWndSyzcsJCaxhr6pPXhilFXcMmwSzgp76SD/8vDHXZvg+qNUF22//uO0Oc97Y0rY5DWe/8wzhvRKpxz2pjWUbVINHXvlIqxbTX1/OndDTz27nq272lk7KAcZk46notP7N9tL6ZC8HbID7d9yIJ1C3h53ctU1lWSmZzJBYMv4JLCSxg/YDxJCaEffVM97Ny0L7Sry6D687DPZdBUt/8GkjMhZxBkF8CxJ0P2oOCrV//9wzq1t46qRboBBX0bSitqeOStdfxlWRn1TQEuGH0M088ZxmlDc7v1efe1O9byYumLvLTuJcpqykhJSOGcARO4JO9kzk7JJ23nVvjoRXj74X1H5TVfHLiirP7BEO//JThu8r4gzy4IvtJzdYeESA+ioA9xd5Zu2M7sRaUsXP0FyYkJXD6ugO+dVciIft30lsLmRjZvXcFLa+fz0pZ3+LSunATgdDL4l9pkzq/aSq+1jwOP71smKX1fYI+8CHIG75vOHgS9jw3eISIiceOoD/rmgPPqqq3MfrOUDz7fQU5GMj88bwTXnDGU/F7dIPDqdkLlmuBrWwlUraWqegOv1W1hQVITy9OCNZ5UV8+tu3dzsaeT16sv5I6CoZOCAd5ymiV7EGT01dG4yFHmqA36PQ1NPLOsjEfeXMfnVXsY0jeDn089kctPjcEF1uZG2L4BKkv2BXoo3JtrvmBjUhIlKcmUpKayMrM37yVDc1Yiw5P78MM+JzF50HkM6j8ueDSe3D2/ZSsisXPUBX3Frnoee3c9f3pvAzv2NHLyoBxum3w8F3X2BVZ3qClvFeZrg9Pb1+OBJrYlJlCSkkJJZg6fZWZT0r8vpYEM6j04zF+CJTCk9xCuHXQelxRewqjcUd36moGIdA9HTdCvKa/h0bdK+cvyTTQ2B7gwdIH11CFRvsDasHtfgG8LnXKpDIV6/U4AdptRkp5JSc4ASnJ7U5I3ljXNNexo3vdtzfz0LEbmjuSbOSMYmTuSkbkjGZ49nLQk3XIoIocmroPe3VmyfjuzF61l4epyUpMSuOLUAr5/ViHD8g/zAmsgALsrYNeW4Gv7+v1OtbBz096ujcD6PoMoyT6GNSPGUZIQoKRxJ5vqq0I96smwREZkj+D8nGCYj8odxYicEeSm5R7p7ouIABEGvZnNAb4KlLv7mDbmG/AAcAmwB/i2uy8PzbsW+PdQ1xnuPi8ahR9MIOC8vGorDy8qpXjjDnIzkrn5/JFcc8YQ8rLaucDqDvW79gX4zi2wazPs2go7Q++7tgRvRww07b9oWjZb+w6jpOBEPksrosSaKWnYzro9W2gKNAHlJNVWMTR7KCflnc7loSP0kbkjGZA5gATTveYi0nkiPaKfC/we+FM78ycDI0Ov04GHgNPNrA/BwcSLCD4tepmZzXf37UdS9MGs2LiD/3phFR98voOhfTOYMW0Ml4/tR3p9BVStgA0tQb6Zxp2bqavZQl3NVupqyqltrqMuwaiz4KvWjLrUTOrSsqlL603tsYXUpZxIXVIqdUkp1CYksrlpF2t2bqCmsRJqK6EWBmQOYGTuSM4ZehEjckcwMmckw7KHkZzYfUaHEpGjR0RB7+6LzGzoQbpMBf7kwZHG3zOzHDMbAEwEXnP3KgAzew2YBDx1JEW3paa+iX+ddw41jQ1kGUwaaZg18fqH9SwobqJ2b4An7A3zJjMwoBfQ62CnShqBSqivJLkxmbSkNNIT00lNSiU/PZ9Lh13KqNxRjMwdyYicEfRK6ZxnSouIHI5onaMfCGwMmy4LtbXXfgAzmw5MBxg8ePAhF5CenEg5NTSnOomJCdRgpCYmk5aaRe+kDI5JySQttTdpqTmkpeWSnp5DWlI6aYlpweBOSictKW3/6dDnlunUxNR9jw4QEekhooRdGFEAAAq3SURBVJVabd224gdpP7DRfTYwG6CoqKjNPgeTmGD85fvF3foZNCIisRCtq4BlwKCw6QJg80HaO4VCXkTkQNEK+vnAP1vQBKDa3bcArwAXmVmumeUCF4XaRESki0R6e+VTBC+s5plZGcE7aZIB3H0WsIDgrZVrCN5e+Z3QvCoz+zmwJLSqu1ouzIqISNeI9K6bqzqY78CN7cybA8w59NJERCQa9E0dEZE4p6AXEYlzCnoRkTinoBcRiXMKehGROKegFxGJcwp6EZE4p6AXEYlzCnoRkTinoBcRiXMKehGROKegFxGJcwp6EZE4p6AXEYlzCnoRkTinoBcRiXMKehGROKegFxGJcxEFvZlNMrNPzWyNmd3axvzfmNmK0OszM9sRNq85bN78aBYvIiId63DMWDNLBB4ELgTKgCVmNt/dP27p4+4/Duv/Q+CUsFXUuvvJ0StZREQORSRH9OOBNe5e6u4NwNPA1IP0vwp4KhrFiYjIkYsk6AcCG8Omy0JtBzCzIUAh8Lew5jQzW2pm75nZtPY2YmbTQ/2WVlRURFCWiIhEIpKgtzbavJ2+VwLPuHtzWNtgdy8CvgX81syGt7Wgu8929yJ3L8rPz4+gLBERiUQkQV8GDAqbLgA2t9P3SlqdtnH3zaH3UuAf7H/+XkREOlkkQb8EGGlmhWaWQjDMD7h7xsyOA3KBd8Pacs0sNfQ5DzgT+Lj1siIi0nk6vOvG3ZvM7AfAK0AiMMfdV5nZXcBSd28J/auAp909/LTOaOBhMwsQ/KPyy/C7dUREpPPZ/rncPRQVFfnSpUtjXYaISI9hZstC10MPoG/GiojEOQW9iEicU9CLiMQ5Bb2ISJxT0IuIxDkFvYhInFPQi4jEOQW9iEicU9CLiMQ5Bb2ISJxT0IuIxDkFvYhInFPQi4jEOQW9iEicU9CLiMQ5Bb2ISJxT0IuIxDkFvYhInIso6M1skpl9amZrzOzWNuZ/28wqzGxF6PX9sHnXmllJ6HVtNIsXEZGOdTg4uJklAg8CFwJlwBIzm9/GIN9/dvcftFq2D3AnUAQ4sCy07PaoVC8iIh2K5Ih+PLDG3UvdvQF4Gpga4fovBl5z96pQuL8GTDq8UkVE5HBEEvQDgY1h02WhttYuN7OVZvaMmQ06xGUxs+lmttTMllZUVERQloiIRCKSoLc22rzV9AvAUHc/CVgIzDuEZYON7rPdvcjdi/Lz8yMoS0REIhFJ0JcBg8KmC4DN4R3cvdLd60OTfwBOjXRZERHpXJEE/RJgpJkVmlkKcCUwP7yDmQ0Im5wCrA59fgW4yMxyzSwXuCjUJiIiXaTDu27cvcnMfkAwoBOBOe6+yszuApa6+3zgJjObAjQBVcC3Q8tWmdnPCf6xALjL3as6YT9ERKQd5t7mKfOYKioq8qVLl8a6DBGRHsPMlrl7UVvz9M1YEZE4p6AXEYlzCnoRkTinoBcRiXMKehGROKegFxGJcwp6EZE4p6AXEYlzCnoRkTinoBcRiXMKehGROKegFxGJcwp6EZE4p6AXEYlzCnoRkTinoBcRiXMKehGROKegFxGJcxEFvZlNMrNPzWyNmd3axvx/NbOPzWylmb1uZkPC5jWb2YrQa37rZUVEpHN1ODi4mSUCDwIXAmXAEjOb7+4fh3X7AChy9z1mdgPw38A3Q/Nq3f3kKNctIiIRiuSIfjywxt1L3b0BeBqYGt7B3f/u7ntCk+8BBdEtU0REDlckQT8Q2Bg2XRZqa8/3gJfCptPMbKmZvWdm0w6jRhEROQIdnroBrI02b7Oj2dVAEXBuWPNgd99sZsOAv5nZh+6+to1lpwPTAQYPHhxBWSIiEolIjujLgEFh0wXA5tadzOwC4A5girvXt7S7++bQeynwD+CUtjbi7rPdvcjdi/Lz8yPeARERObhIgn4JMNLMCs0sBbgS2O/uGTM7BXiYYMiXh7Xnmllq6HMecCYQfhFXREQ6WYenbty9ycx+ALwCJAJz3H2Vmd0FLHX3+cB9QBbw/5kZwOfuPgUYDTxsZgGCf1R+2epuHRER6WTm3ubp9pgqKirypUuXxroMEZEew8yWuXtRW/P0zVgRkTinoBcRiXMKehGROKegFxGJcwp6EZE4p6AXEYlzCnoRkTinoBcRiXMKehGROKegFxGJcwp6EZE4p6AXEYlzCnoRkTinoBcRiXMKehGROKegFxGJcwp6EZE4p6AXEYlzCnoRkTgXUdCb2SQz+9TM1pjZrW3MTzWzP4fmLzazoWHzbgu1f2pmF0evdBERiUSHQW9micCDwGTgBOAqMzuhVbfvAdvdfQTwG+De0LInAFcCJwKTgJmh9YmISBeJ5Ih+PLDG3UvdvQF4Gpjaqs9UYF7o8zPA+WZmofan3b3e3dcBa0LrExGRLpIUQZ+BwMaw6TLg9Pb6uHuTmVUDfUPt77VadmBbGzGz6cD00GSNmX0aQW1tyQO2HeayPZX2Of4dbfsL2udDNaS9GZEEvbXR5hH2iWTZYKP7bGB2BPUclJktdfeiI11PT6J9jn9H2/6C9jmaIjl1UwYMCpsuADa318fMkoBsoCrCZUVEpBNFEvRLgJFmVmhmKQQvrs5v1Wc+cG3o8xXA39zdQ+1Xhu7KKQRGAu9Hp3QREYlEh6duQufcfwC8AiQCc9x9lZndBSx19/nAo8BjZraG4JH8laFlV5nZ/wAfA03Aje7e3En70uKIT//0QNrn+He07S9on6PGggfeIiISr/TNWBGROKegFxGJc3ET9B09piHemNkgM/u7ma02s1VmdnOsa+oqZpZoZh+Y2f/GupauYGY5ZvaMmX0S+nmfEeuaOpuZ/Tj0e/2RmT1lZmmxrinazGyOmZWb2UdhbX3M7DUzKwm950ZjW3ER9BE+piHeNAH/j7uPBiYANx4F+9ziZmB1rIvoQg8AL7v78cBY4nzfzWwgcBNQ5O5jCN4EcmVsq+oUcwk+GibcrcDr7j4SeD00fcTiIuiJ7DENccXdt7j78tDnXQT/52/zW8fxxMwKgEuBR2JdS1cws97AOQTvbMPdG9x9R2yr6hJJQHroezkZxOH3b9x9EcG7FMOFP05mHjAtGtuKl6Bv6zENcR96LUJPCz0FWBzbSrrEb4GfAIFYF9JFhgEVwB9Dp6seMbPMWBfVmdx9E/Ar4HNgC1Dt7q/Gtqouc4y7b4HgwRzQLxorjZegj/hRC/HGzLKAvwA/cvedsa6nM5nZV4Fyd18W61q6UBIwDnjI3U8BdhOlf853V6Hz0lOBQuBYINPMro5tVT1bvAT9UfmoBTNLJhjyT7j7s7GupwucCUwxs/UET899xcwej21Jna4MKHP3ln+tPUMw+OPZBcA6d69w90bgWeDLMa6pq3xhZgMAQu/l0VhpvAR9JI9piCuhx0A/Cqx29/tjXU9XcPfb3L3A3YcS/Bn/zd3j+kjP3bcCG83suFDT+QS/aR7PPgcmmFlG6Pf8fOL8AnSY8MfJXAv8NRorjeTpld1ee49piHFZne1M4BrgQzNbEWq73d0XxLAm6Rw/BJ4IHcSUAt+JcT2dyt0Xm9kzwHKCd5d9QBw+DsHMngImAnlmVgbcCfwS+B8z+x7BP3hfj8q29AgEEZH4Fi+nbkREpB0KehGROKegFxGJcwp6EZE4p6AXEYlzCnoRkTinoBcRiXP/P+/GbmMewUz+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(sol.t, h1,\n",
    "         sol.t, h2, \n",
    "         sol.t, h3)\n",
    "plt.ylim(0, 2)\n",
    "plt.legend(['$h_1$', '$h_2$', '$h_3$'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1519bf0e10>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de3Scd33n8fdXM9KMNJKty4xs2bJkObFl2bk5VmIngZALhiRAHLqBJoFu6bLJCYfUXQptKCekwLKHtLvApmxOOWlg4ewCgVKauG3aAC3ZXkgcy07S4GuMr4ovutuy7qP57R/PaDSWxpacSHpmRp/XOXM0z/P8ZvQdsD/5+Tu/53nMOYeIiOS+Ar8LEBGRmaFAFxHJEwp0EZE8oUAXEckTCnQRkTwR9OsXR6NRt3z5cr9+vYhITtqxY0eHcy6W6Zhvgb58+XJaWlr8+vUiIjnJzI6c75haLiIieUKBLiKSJxToIiJ5QoEuIpInFOgiInlCgS4ikicU6CIiecK3degiIrnCOUc84RiOJ7zHqPdzaMK293z0vMdGks9vbVrElcvKZ7xOBbqIZJ346LlBODQpNBOMxBMMpe/LMGYoPjlopw5k773Pef1ogpm8dUT1grACXURm32jC0Tcc5+xgnN7BOGeH4gzFJ4fhlME4tj1F6I6/fjS1nZjB8CwKFFAUTD4mPC8MFhAKFBApClJRMuF4ckwo02vTtsePByYcs3P3jf3OgGFmM/cB0yjQRfJEYiyIh7wwPjM4/rx3cISzQ+MBnb7dO2Fc3/DoW64hWGDnDc+x4AsXFrAgHEweC2QIxvHXFY4dO18oTwzdDL93tsIzGynQRXzmnKN/eDQZrCOTQvZMMnzPpsI4Tu9QnLODI+eMOzscn1ZbIFIUoCxcSGk4SGkoSFk4yJLyMKWhIKWhQsrC3j7vWCGRUIBQMJAxdAvTwzhQQEHB/AnPbKRAF3mLnHMMjIxmnA33DqW3LM4/G+4ditM3FJ9Wi6GkKJAK4NJwIWWhIIsWJIM4HKQsGcDpQe0Fs7evLBwkUhQkoNDNWwp0mZecc/QOxenoHUprQ4y3IsZmw+NBPZJx3HSCuLgwkBa4XvhGS0u88A2dOxueOK40FKQs5M2SgwGtMpYLU6BLXkkkHF39w7SdGaKtd5C23iHae4doO5P2vNc7NjiSuOB7hQsLUi2IseCti5RQGg6yIBnGmWbD6eNLQ0EFscwZBbrkhOF4go6zyTBOhrMX0IPJ8PbCuuPsEPEM0+aycJDqshDVZWHW1ZVTXRYilnykh/PY80jI+9JOJJco0MVXfUPxtJnzeDi39Q4mZ9be8+7+kUmvNYOqSBGxsjDVZSFWLy6jeoEX2rGyUCrAY2UhiosCPnw6kbmlQJcZ55yjp3/k3GBOC+f0NkimJXKFASNWGiK2IExdVQnNyyuoLgsnw3o8pKOlRWpniKRRoMu0OedoPzvEqdMTgnlC26O9d4jh0cn96UhRgOoFYWKlIdYsWcBNjTEvqMtCqZl1dVmI8pLCebV2WGSmKNBlktP9Ixzq7ONQx1kOtfdxqLM/9TzTjLqipDA1g14RiyTbHcmgLgtRvcB7Hgnpj5vIbNLfsHlqYHiUw519HOqY/OjqG06NKzCorSihIRqhub6S5VUl1JQXp4I6VhrSl4ciWUKBnsdGRhMc6+qfFNiHO/o4fnrwnLGLFoRoiEZ479pFNEQjNERLaYhGWFZZTCioLxRFcsG0At3MbgMeBwLAU865xyYcrwO+C5Qnx3zWOffcDNcqGSQSjhNnBjnc0cfBjj6vRdJxlsOd/Rzt6mc0bQnfwuJCGqIRNq6o8kI7FmF5VYTl0QilaoeI5Lwp/xabWQB4AtgEtALbzWyrc2532rBHgB855/7czNYAzwHLZ6Heeck5R1ffMIfGQjs5yx6bcQ/Fx7+ALC4MsDwaYU3NAt53eQ3LoxEaohFWRCNURIp8/BQiMtumMy27FjjgnDsIYGZPA5uB9EB3wILk84XA8Zkscj5p7x1ix5Fu9p3s9b6ITIZ472A8NSZYYNRVldBQFeEdl0ZpiEWSbZIIi8rCukCSyDw1nUBfChxL224FNkwY8wXgp2b2u0AEeHemNzKzB4AHAOrq6i621rwzmnC80dZLy+Fudh7ppuVIN0e7+gHvpJklC4tpiEa466qlqcBuiEaorSjW+msRmWQ6gZ5pujfx3Op7ge84575qZtcB/8fMLnPOnbMY2Tn3JPAkQHNz8wxewj43nB2K8+rRHnYc6ablSBevHu2hd8ibeUdLi1hfX8FHN9axvr6StUsWEC7Ul5EiMn3TCfRWYFnadi2TWyofB24DcM69aGZhIAq0zUSRucg5R2v3ADuPdtNyuJsdR7rZe/IMCefNvhsXlfGBq5bQXF/B+voK6ipLdDKNiLwt0wn07cBKM2sA3gTuAe6bMOYocCvwHTNrAsJA+0wWmu2G4wl2nzhDy+Eudh71AvzUmSHAO0NyXV0FD92ykvX1FayrK2dBuNDnikUk30wZ6M65uJk9BDyPtyTx2865XWb2JaDFObcV+DTwF2b2Kbx2zMecm8lbqmaf7r5hb/Z9xAvv1471pFabLC0vZuOKKtYnZ9+Ni8rU8xaRWWd+5W5zc7NraWnx5XdfLOccv27vY8eRrmT/u5uD7X2At+Jk7dKFrK+roHl5BVfXVbB4YdjnikUkX5nZDudcc6ZjOpvkAvaf6uX7246y9bXjqdPhy0sKWV9Xwd3ra1lfV8EVteW6NKuIZAUF+gSDI6P8w69O8v1tR3n5cBdFgQLes3YRN66McXV9BZfEIvryUkSykgI96VBHHz94+Sh/2XKM7v4R6qtK+Nwdq7l7/TIqdYaliOSAeR3oI6MJfr77FN/bdpR/PdBBoMB4z5pFfGRDPddfUqUzLkUkp8zLQG/t7ufpl4/xw5ZjtPcOsWRhmE9vWsWHr1nGogX6QlNEctO8CfTRhOOFfW18b9tRfrHPO9/plsZqPrKxjnetqiag2biI5Li8D/S2M4P8cPsxfvDyUY6fHiRWFuKhmy/lN69ZRm1Fid/liYjMmLwM9ETC8W+/7uB7Lx3lZ3tOMZpwvOPSKJ9//xrevWYRhTrJR0TyUF4FeufZIX68o5Xvv3yUI539VJQU8p/f0cC919axPBrxuzwRkVmV84HunGP74W6+t+0If//6SYZHE1y7vJLf37SK965drCsWisi8kbOBnkg4/mpnK0/+80HeaDtLWTjIfRvquG9DHasWlfldnojInMvJQN95tJsvbt3Fa62nuWzpAv707iv4wBVLdAq+iMxrORfoT/3LQb78d3uoLgvx9d+8kruuWqpT8UVEyMFAv/6SKA++6xIeuuVS3aleRCRNziXimiULWLNkwdQDRUTmGS3IFhHJEwp0EZE8oUAXEckTCnQRkTyhQBcRyRMKdBGRPKFAFxHJEwp0EZE8oUAXEckTCnQRkTyhQBcRyRMKdBGRPKFAFxHJEwp0EZE8oUAXEckTCnQRkTyhQBcRyRMKdBGRPKFAFxHJEwp0EZE8Ma1AN7PbzGyfmR0ws89mOP51M3s1+dhvZj0zX6qIiFxIcKoBZhYAngA2Aa3AdjPb6pzbPTbGOfeptPG/C6ybhVpFROQCpjNDvxY44Jw76JwbBp4GNl9g/L3AD2aiOBERmb7pBPpS4Fjadmty3yRmVg80AP90nuMPmFmLmbW0t7dfbK0iInIB0wl0y7DPnWfsPcCPnXOjmQ465550zjU755pjsdh0axQRkWmYTqC3AsvStmuB4+cZew9qt4iI+GI6gb4dWGlmDWZWhBfaWycOMrNGoAJ4cWZLFBGR6Zgy0J1zceAh4HlgD/Aj59wuM/uSmd2ZNvRe4Gnn3PnaMSIiMoumXLYI4Jx7Dnhuwr5HJ2x/YebKEhGRi6UzRUVE8oQCXUQkTyjQRUTyhAJdRCRPKNBFRPLEtFa5zJWRkRFaW1sZHBz0u5SLFg6Hqa2tpbCw0O9SRGSeyqpAb21tpaysjOXLl2OW6YoD2ck5R2dnJ62trTQ0NPhdjojMU1nVchkcHKSqqiqnwhzAzKiqqsrJf1mISP7IqkAHci7Mx+Rq3SKSP7Iu0EVE5K1RoIuI5AkFuohInlCgZ3DzzTfzs5/9DIBHHnmELVu2+FyRiMjUsmrZYrov/s0udh8/M6PvuWbJAv74A2un/t1f/CKPPvoobW1tvPLKK2zdOuny7yIiWUcz9AxuvPFGnHN87Wtf4+mnnyYQCHDw4EE+/vGPc/fdd/tdnohIRlk7Q5/OTHq2vP7665w4cYJoNEpZWRkAK1as4Fvf+pYCXUSylmboE5w4cYKPfOQjPPvss0QiEZ5//nm/SxIRmRYFepr+/n5+4zd+g69+9as0NTXx+c9/ni984Qt+lyUiMi0K9DQlJSW8+OKLbNq0CfB66S++6N3zurOzkwcffJBXXnmFr3zlK36WKSKSUdb20LNNVVUV3/zmN/0uQ0TkvDRDFxHJEwp0EZE8oUAXEckTCnQRkTyhQBcRyRMKdBGRPKFAFxHJEwp0EZE8oUAXEckTCvQMdIMLEclF2Xvq/99/Fk6+PrPvufhyuP2xKYfpBhcikos0Q88g0w0unnnmGe6//342b97MT3/6U79LFBGZJHtn6NOYSc+WTDe4uOuuu7jrrrvo7u7mM5/5DO95z3t8q09EJBPN0CeY6gYXX/7yl/nkJz/pU3UiIuc3rUA3s9vMbJ+ZHTCzz55nzIfNbLeZ7TKz789smXPjQje4cM7x8MMPc/vtt3P11Vf7W6iISAZTtlzMLAA8AWwCWoHtZrbVObc7bcxK4I+AG5xz3WZWPVsFz6axG1yMSb/BxTe+8Q1+/vOfc/r0aQ4cOMCDDz7oV5kiIhlNp4d+LXDAOXcQwMyeBjYDu9PG3A884ZzrBnDOtc10oX7bsmWLli+KSFabTstlKXAsbbs1uS/dKmCVmf2bmb1kZrdleiMze8DMWsyspb29/a1VLCIiGU0n0C3DPjdhOwisBG4C7gWeMrPySS9y7knnXLNzrjkWi11srSIicgHTCfRWYFnadi1wPMOYZ51zI865Q8A+vIAXEZE5Mp1A3w6sNLMGMysC7gEmnjr5DHAzgJlF8VowB2eyUBERubApA905FwceAp4H9gA/cs7tMrMvmdmdyWHPA51mthv4BfAHzrnO2SpaREQmm9aZos6554DnJux7NO25A34/+RARkTHOQXwQRgZgpB9GBiFSBcUVM/6rsvfUfxGR2TQ6kgzYgXMf8QnbI/3JQJ7O2PR9yfCOD0z+3e//OjT/pxn/SAp0EckeicTkkJwYkGPBmgrZ9H0ZAvacQE6bKbvRt1CgQWEJFBYnf4bHnxdFIBI7d18wnDY++QgWw9LZOdtcgS4iF+YcjA5nCMgJITvlLHbi2AwhPTr01moMhM4NzVSYFkPpomTITgzYtOfB4gkhnfb69H2BIrBMK7mzQ9YG+p+8/Cfs7do7o++5unI1D1/78JTjbr75Zj73uc+xadMmHnnkEc6cOcOf/dmfzWgtIr4a7oe+dujrSP5sh762Cdsd4z/fymzWCs6dnaaHZqgsGbTFaSGbKWAzhfTEGW8YCgIz/79RDsraQPeTbnAhOScxCv1daWE88ZEM57PJ0B7py/w+RWUQiXqtg/J6WLoeSqq8dsLEIJ1qxhsozOrZbD7K2kCfzkx6tqTf4OKFF14gEAiwZ88eHn/8cTo6Orj11lv5xCc+4Vt9Mg84B8N90wvnvnbo72TyCdyABbxwjsS8oF7W4D0vjZ27PxKDkigUlcz5R5WZk7WB7qdMN7hoamrim9/8JolEgvvvv9/nCiUnjcahP1NLI/k4O2F/ptURAKGF4yFcdQnUXzc5nMce4XIo0G0P5gsF+gTpN7jYsmULzz//PO9973sB2Lp1K4899hgPPfSQz1VKVnAOhs5MDua+juTseUJoD3Rlfp+CwnPDONY4OZgjUYhUez+Dobn9nJIzFOhpMt3g4uGHH04F+p133smdd97J+973Pu677z6fq5VZER+ePIvOFM5jz8+3KiNcnmxtVEP1aoi8M8MMOhnQ4YXqNcuMUKCnudANLl544QV+8pOfMDQ0xB133OFXifJ2xYeg8wC07fEenW+ktTraYPB05tcFQl44R6Lez0VrM8yik4+SKggWze3nEkGBPm033XQTN910k99lyHSNxqHrILTthva93s+2vV6Yjy3BswBUNkBZDSy+LHM4j4V2qEyzaMl6CnTJbYkE9Bwen3G37fECvGO/dzIMAOYFd6wJmj4A1U3eo+pS9aMlryjQJTc4B6dbk4GdHt77zl0NsrDO61lfeqsX4NVNEF2l5XgyLyjQJbs4B2dPjbdIUi2TvTDcOz6urAZiq70LHFWvhuo13uqQUJl/tYv4LOsC3TmH5WCv0ruCsFyUvs7JPe623TDYMz6mJOrNsq+61wvw6jVegM/CpUdFcl1WBXo4HKazs5OqqqqcCnXnHJ2dnYTDYb9LyU4DPcnQ3nNuy6Qv7Ubh4YVei2TtB8d73LEm74xGEZmWrAr02tpaWltbaW9vn3pwlgmHw9TW1vpdhr+Gzno97fQed9se6E27BW1RqdcaWfXe8R53dZPXQsmh/4iLZKOsCvTCwkIaGhr8LkOmMjLgrSI5p8e9G3qOjo8Jhr0vIxtuTOtxr4aFy3QqusgsyapAlyw10AOH/wV+/QvvZ+cBcAnvWEEhRFdC7TWw7j+Oz7grluuSpiJzTIEuk8WHofVlOPiCF+LHd3oBXhiB+uvH+9yxJu/iUIFCvysWERToAt5Swbbd4wF+5N+8O8hYwLse9js/Aytu8mbhOqVdJGsp0OerM8fHA/zgC951TMA7e/Kqj3gB3vBOb/WJiOQEBfp8MdQLh/91PMQ79nn7S6JeeI89ypf5VKCIvF0K9Hw1OgJv7oSDyRl463ZIxL3VJ/XXw7qPwiU3Q/VarToRyRMK9HzhHHS8MR7gh/4leaq8wZKr4Pot3gx82Qbv3o8ikncU6LnsbBsc/H/jIX7mTW9/xXK4/O5kH/xGKKn0r0YRmTMK9Fwy3A9Hfjke4Kd+5e0Pl8OKd8GKP/BCvFInZ4nMRwr0bJYYheOvjgf4sW3eNb4DRVC3EW79Yy/Aa67USTwiokDPOn2dsO/v4I2fwqF/Hr8l2uLLYcODXoDXXafre4vIJAr0bNDXAXv+BnY/432Z6UZhQa13d50VN0PDu3TVQRGZkgLdL2fbYM9W2P2stz7cJaByBdzwe7D2Llh8ha4+KCIXRYE+l3pPejPxXc94p9fjoGolvPPTsGYzLLpMIS4ib5kCfbadOQ67kzPxoy8CzruM7Lv+ENbc5V3kSiEuIjNAgT4bTrcmQ/wZb2UKeNcDv+mPvJl49Wp/6xORvKRAnyk9R71Z+O5nvdPswWuh3PyIF+KxVf7WJyJ5b1qBbma3AY8DAeAp59xjE45/DPjvQPJURf6Xc+6pGawzO3Uf9gJ81zPeNcPB+zLz1kehaTNEL/W1PBGZX6YMdDMLAE8Am4BWYLuZbXXO7Z4w9IfOuYdmocbs0nVwPMRPvOrtq7kK3v0FbyZeucLP6kRkHpvODP1a4IBz7iCAmT0NbAYmBnp+cg46f+31w3c/Aydf9/YvXQ+bvuSFeMVyX0sUEYHpBfpS4FjadiuwIcO4/2BmNwL7gU85545NHGBmDwAPANTV1V18tXMhkfDuWn/0RTjyovdz7KJXtdfAe/4brLkTyrO0fhGZt6YT6JnW1LkJ238D/MA5N2RmDwLfBW6Z9CLnngSeBGhubp74Hv6ID8PxV+DoL+HoS95jsMc7VroY6q+D+hug8XZYWOtvrSIiFzCdQG8F0m9jUwscTx/gnOtM2/wL4E/efmmzZPCMdwPksdn3mzsgPugdq1rpnW5ff713vZSK5VojLiI5YzqBvh1YaWYNeKtY7gHuSx9gZjXOuRPJzTuBPTNa5dvRe8oL7qMvepeePfUr7zR7C0DNFdD8ce/KhXXX6XopIpLTpgx051zczB4Cnsdbtvht59wuM/sS0OKc2wpsMbM7gTjQBXxsFmu+ULHeKpRU//uX3jZAsBhqm+HGP/DCu/YaCJX6UqaIyGww5/xpZTc3N7uWlpa3/gajIzDSD12Hkr3vZA/87CnveHGFF9x113ktlMVXQLBoZooXEfGJme1wzjVnOpZ7Z4r+8hvws0e9tkm6hXXeZWbrr4O66yG6Sjc/FpF5JfcCfel6eMfve3evLyyGssVeD1wrUERknsu9QK+/3nuIiMg5ci/QRUR8Fk/EGYgPpB79I/3ez3j/OdupfSPjxwbiA3xo1Ye4YekNM16XAl1E8pJzjuHE8KQwnTJ8J44ZmXxsODF8UbUUB4vPefSO9M7KZ1agi4ivEi7BYHxw0kw208w2474JAZ1+fNSNTruOgAUoCZZQHCympLAkFb4LwwupCdaktseOjY3NtC99OxwMU2Bzs0BDgS4i0zKSGLlgeE46NkXYpj8uRigQyhiei0sWewFbmDlYpwrfwoJCLMfPDFegi8xjQ6NDtPe30z7QTlt/Gx0DHbT1t6X2tfe30zHYQd9IH/FEfNrva9g5s9yx8CwtKiVWEjt/wBZeOHyLg8UECgKz+L9IblOgi+Sh4dHhVCCPhXV6SI/tOzN8ZtJrgwVBqouriZZEaVjYQPPiZkoLSy84u514LBQI5fxsNxcp0EVyyMjoiDeLHvAC+pxZ9cB4YPcM9Ux6bdCCREuiVBdXU7+gnvWL1lNdUk2sOEZ1STXR4ijVJdWUh8oVxjlKgS6SBUYSI3QOdI4HczKs00O6vb+d7qHuSa8NWIBocZRYcYza0lqurr6aWHGMWEksFdaxkhjlofI5+3JO/KFAF5lFY0Gd3ptuG5jcq+4a7Jr02gIrIBqOEiuJsaR0CVfGriRWEqO6uDoV1rGSGJXhSgW1AAp0kbcknojTNdiVcSad3gbpGuzCTbgfTIEVUBWuIlocZXFkMZfHLk/1rMfCurqkmopQhb4AlIuiQBdJM5oYpWuwy5tF94/3qie2QboGu0hMuECcYVSGK1MtjjVVa85te6TNqIMF+qsnM09/qmTeOTt8lv3d+9nbtZdf9/yatv62VIB3DHZMCmogFdTR4ihNVU2TvkiMFceoKq5SUIuv9KdP8pZzjpN9J9nXvY+9XXvZ1+X9bD3bmhpTVlRGTaSGWEmMxorG8YBO61VXhasoDBT6+ElEpkeBLnlhZHSEg6cPTgrv9HXW9Qvqaapq4oMrP8jqytU0VjRSXVKtJXqSNxToknNOD51mf/f+VGjv697HgZ4DqTMZQ4EQqypWsal+E6srV7O6cjUrK1YSKYz4XLnI7FKgS9ZyzvHm2TfZ17XvnJn38b7jqTGV4UqaKpu4bs11rK7wwrtuQZ162TIv6U+9ZIXh0WEO9Bw4J7z3d+1PXWbUMJYvXM6VsSv5UOOHUjPvaHHU58pFsocCXeZcz2APe7u92fa+rn3s7d7LoZ5DxJ3XMikOFrOqYhV3rLiDxspGGisaubT8UkoKS3yuXCS7KdBl1iRcgtbe1lSfe6znfar/VGpMdXE1jZWN3FR7Uyq8l5Ut0wk1Im+BAl1mxGB8kAM9B1J97rEA74/3A971Rsau3Le6YrUX3pWNVIYrfa5cJH8o0OWidQ50plolY22TQ2cOpU7IiRRGaKxoZPOlm73lgZVeyyQUCPlcuUh+U6DLlE72neSlEy+x7cQ2Xj75Mm39baljNZEaGisaeXf9u1PhvbR0qS4WJeIDBbpM0jPYw8snX2bbiW1sO7mNI2eOAN4SwWsXX8vl0ctT4b0wtNDnakVkjAJdGIgPsPPUTrad2MZLJ15ib9deHI6SYAnNi5v58KoPs6FmAysrVmrmLZLFFOjz0EhihF0du1JtlFfbXyWeiBMsCHJl7Eo+cdUnuK7mOtZG11JYoGuYiOQKBfo84JzjjZ43UjPwlpMt9Mf7MYzVlav5rabfYkPNBtZVr9Nab5EcpkDPU629rV4PPNkHH7sjTv2Cet6/4v1sqNnAtYuvpTxc7nOlIjJTFOh5omuwi5dPvJxqo4xdIjZWHOP6JdezoWYDG2s2sjiy2OdKRWS2KNBzVP9IPy2nWlIBvr97PwBlhWU0L27mo2s+ysaajaxYuEKXhxWZJxToOWJkdITX2l9j20mvjfJ6++vEXZyigiLWVa9jy7otbKzZSFNVk640KDJP6W9+lkq4BPu69nlfZJ58iZ2ndjIQH6DAClhbtZaPXfYxNtRs4KrYVYSDYb/LFZEsMK1AN7PbgMeBAPCUc+6x84y7G/hL4BrnXMuMVTkPjCZG2d+9n51tO9lxagfbT26nZ6gHgBULV3DXpXexoWYD1yy+hgVFC3yuVkSy0ZSBbmYB4AlgE9AKbDezrc653RPGlQFbgG2zUWi+GRkdYVfnLlpOtbDj1A5ebXuVsyNnAVhaupQba29kY81GNtRsoLqk2udqRSQXTGeGfi1wwDl3EMDMngY2A7snjPuvwJ8Cn5nRCvNE/0g/r7W/xo5TO9jZtpN/b/93hkaHALhk4SXc3nA76xetZ/2i9VqJIiJvyXQCfSlwLG27FdiQPsDM1gHLnHN/a2bnDXQzewB4AKCuru7iq80hp4dOs/PUzlQLZU/nHuIuToEV0FTZxIcbP8z6Reu5uvpqKsIVfpcrInlgOoGeac2bSx00KwC+Dnxsqjdyzj0JPAnQ3NzsphieU9r629h5ameqhXKg5wAARQVFXBa9jN+57HdoXtTMldVX6mbFIjIrphPorcCytO1a4HjadhlwGfBCcr3zYmCrmd2Zr1+MOudo7W1NhffOtp0c6/X+EVMSLGFd9TruaLiDqxddzWXRy3QdcBGZE9MJ9O3ASjNrAN4E7gHuGzvonDsNpO7Ua2YvAJ/JpzBPuAQHeg544X3Ka6G0D7QDUB4qZ/2i9dzTeA/rF6+nsaJR68BFxBdTJo9zLm5mDwHP4y1b/LZzbpeZfQlocc5tne0iZ8NIYoSewR66h7rpHvQeXYNdqe2uwS56hnroHuzmVN+p1N3nF5Us4prF16S+wNSZmCKSLcw5f1rZzc3NrqVl5ibxg/FBL4iHus4N6MFueoZ6Us+7h7z9vcO9GYApJeMAAAQ6SURBVN/HMBaGFlIRrqAiVEFluJKq4iquiF3B+kXrWRJZogAXEd+Y2Q7nXHOmYznXG/jrN/6a7+z6Tmp7MD5I91A3A/GBjOODFqQiXEF5uJzKUCVNlU1eWIcrqAxVpp5XhLyf5aFy3XFeRHJSzgX6wtBCLim/BEsuvikKFHnhHK5MhXJleDyoywrLNKMWkXkh5wL9lrpbuKXuFr/LEBHJOrpBpIhInlCgi4jkCQW6iEieUKCLiOQJBbqISJ5QoIuI5AkFuohInlCgi4jkCd+u5WJm7cCRt/jyKNAxg+XkAn3m+UGfeX54O5+53jkXy3TAt0B/O8ys5XwXp8lX+szzgz7z/DBbn1ktFxGRPKFAFxHJE7ka6E/6XYAP9JnnB33m+WFWPnNO9tBFRGSyXJ2hi4jIBAp0EZE8kXOBbma3mdk+MztgZp/1u57ZZmbLzOwXZrbHzHaZ2e/5XdNcMLOAmb1iZn/rdy1zwczKzezHZrY3+f/1dX7XNNvM7FPJP9O/MrMfmFnY75pmmpl928zazOxXafsqzexnZvZG8mfFTP2+nAp0MwsATwC3A2uAe81sjb9Vzbo48GnnXBOwEfjkPPjMAL8H7PG7iDn0OPAPzrnVwJXk+Wc3s6XAFqDZOXcZEADu8beqWfEd4LYJ+z4L/KNzbiXwj8ntGZFTgQ5cCxxwzh10zg0DTwObfa5pVjnnTjjndiaf9+L9RV/qb1Wzy8xqgfcBT/ldy1wwswXAjcC3AJxzw865Hn+rmhNBoNjMgkAJcNznemacc+6fga4JuzcD300+/y5w10z9vlwL9KXAsbTtVvI83NKZ2XJgHbDN30pm3f8E/hBI+F3IHFkBtAP/O9lmesrMIn4XNZucc28C/wM4CpwATjvnfupvVXNmkXPuBHgTNqB6pt441wLdMuybF+suzawU+Cvgvzjnzvhdz2wxs/cDbc65HX7XMoeCwNXAnzvn1gF9zOA/w7NRsm+8GWgAlgARM/uov1XlvlwL9FZgWdp2LXn4z7SJzKwQL8y/55z7id/1zLIbgDvN7DBeS+0WM/u//pY061qBVufc2L+8fowX8Pns3cAh51y7c24E+Alwvc81zZVTZlYDkPzZNlNvnGuBvh1YaWYNZlaE9yXKVp9rmlVmZni91T3Oua/5Xc9sc879kXOu1jm3HO//339yzuX1zM05dxI4ZmaNyV23Art9LGkuHAU2mllJ8s/4reT5F8FptgK/nXz+28CzM/XGwZl6o7ngnIub2UPA83jfin/bObfL57Jm2w3AbwGvm9mryX2fc84952NNMvN+F/hecqJyEPgdn+uZVc65bWb2Y2An3kquV8jDSwCY2Q+Am4CombUCfww8BvzIzD6O9x+2D83Y79Op/yIi+SHXWi4iInIeCnQRkTyhQBcRyRMKdBGRPKFAFxHJEwp0EZE8oUAXEckT/x/H4s7ASrKnrwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(sol.t, x1,\n",
    "         sol.t, x2, \n",
    "         sol.t, x3)\n",
    "plt.legend(['$x_1$', '$x_2$', '$x_3$'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}